Copyright © 2024 Daniel Zimmermann.

You may use, distribute and modify this code under the terms of the MIT license.

You should have received a copy of the MIT license with this file. If not, please visit: https://opensource.org/license/mit

Setup

source("helpers.R")
## Lade nötiges Paket: S4Vectors
## Lade nötiges Paket: stats4
## Lade nötiges Paket: BiocGenerics
## 
## Attache Paket: 'BiocGenerics'
## Die folgenden Objekte sind maskiert von 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## Die folgenden Objekte sind maskiert von 'package:base':
## 
##     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
##     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
##     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
##     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
##     Position, rank, rbind, Reduce, rownames, sapply, saveRDS, setdiff,
##     table, tapply, union, unique, unsplit, which.max, which.min
## 
## Attache Paket: 'S4Vectors'
## Das folgende Objekt ist maskiert 'package:utils':
## 
##     findMatches
## Die folgenden Objekte sind maskiert von 'package:base':
## 
##     expand.grid, I, unname
## Lade nötiges Paket: IRanges
## 
## Attache Paket: 'IRanges'
## Das folgende Objekt ist maskiert 'package:grDevices':
## 
##     windows
## Lade nötiges Paket: GenomicRanges
## Lade nötiges Paket: GenomeInfoDb
## Lade nötiges Paket: SummarizedExperiment
## Lade nötiges Paket: MatrixGenerics
## Lade nötiges Paket: matrixStats
## 
## Attache Paket: 'MatrixGenerics'
## Die folgenden Objekte sind maskiert von 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars
## Lade nötiges Paket: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## 
## Attache Paket: 'Biobase'
## Das folgende Objekt ist maskiert 'package:MatrixGenerics':
## 
##     rowMedians
## Die folgenden Objekte sind maskiert von 'package:matrixStats':
## 
##     anyMissing, rowMedians
## 
## Attache Paket: 'dplyr'
## Das folgende Objekt ist maskiert 'package:Biobase':
## 
##     combine
## Das folgende Objekt ist maskiert 'package:matrixStats':
## 
##     count
## Die folgenden Objekte sind maskiert von 'package:GenomicRanges':
## 
##     intersect, setdiff, union
## Das folgende Objekt ist maskiert 'package:GenomeInfoDb':
## 
##     intersect
## Die folgenden Objekte sind maskiert von 'package:IRanges':
## 
##     collapse, desc, intersect, setdiff, slice, union
## Die folgenden Objekte sind maskiert von 'package:S4Vectors':
## 
##     first, intersect, rename, setdiff, setequal, union
## Die folgenden Objekte sind maskiert von 'package:BiocGenerics':
## 
##     combine, intersect, setdiff, union
## Die folgenden Objekte sind maskiert von 'package:stats':
## 
##     filter, lag
## Die folgenden Objekte sind maskiert von 'package:base':
## 
##     intersect, setdiff, setequal, union
## 
## Attache Paket: 'magrittr'
## Das folgende Objekt ist maskiert 'package:GenomicRanges':
## 
##     subtract
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ readr     2.1.5
## ✔ ggplot2   3.5.1     ✔ stringr   1.5.1
## ✔ lubridate 1.9.4     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ lubridate::%within%() masks IRanges::%within%()
## ✖ dplyr::collapse()     masks IRanges::collapse()
## ✖ dplyr::combine()      masks Biobase::combine(), BiocGenerics::combine()
## ✖ dplyr::count()        masks matrixStats::count()
## ✖ dplyr::desc()         masks IRanges::desc()
## ✖ tidyr::expand()       masks S4Vectors::expand()
## ✖ tidyr::extract()      masks magrittr::extract()
## ✖ dplyr::filter()       masks stats::filter()
## ✖ dplyr::first()        masks S4Vectors::first()
## ✖ dplyr::lag()          masks stats::lag()
## ✖ ggplot2::Position()   masks BiocGenerics::Position(), base::Position()
## ✖ purrr::reduce()       masks GenomicRanges::reduce(), IRanges::reduce()
## ✖ dplyr::rename()       masks S4Vectors::rename()
## ✖ lubridate::second()   masks S4Vectors::second()
## ✖ lubridate::second<-() masks S4Vectors::second<-()
## ✖ purrr::set_names()    masks magrittr::set_names()
## ✖ dplyr::slice()        masks IRanges::slice()
## ✖ magrittr::subtract()  masks GenomicRanges::subtract()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Lade nötiges Paket: AnnotationDbi
## 
## 
## Attache Paket: 'AnnotationDbi'
## 
## 
## Das folgende Objekt ist maskiert 'package:dplyr':
## 
##     select
## 
## 
## 
## 
## Lade nötiges Paket: data.table
## 
## 
## Attache Paket: 'data.table'
## 
## 
## Die folgenden Objekte sind maskiert von 'package:lubridate':
## 
##     hour, isoweek, mday, minute, month, quarter, second, wday, week,
##     yday, year
## 
## 
## Das folgende Objekt ist maskiert 'package:purrr':
## 
##     transpose
## 
## 
## Die folgenden Objekte sind maskiert von 'package:dplyr':
## 
##     between, first, last
## 
## 
## Das folgende Objekt ist maskiert 'package:SummarizedExperiment':
## 
##     shift
## 
## 
## Das folgende Objekt ist maskiert 'package:GenomicRanges':
## 
##     shift
## 
## 
## Das folgende Objekt ist maskiert 'package:IRanges':
## 
##     shift
## 
## 
## Die folgenden Objekte sind maskiert von 'package:S4Vectors':
## 
##     first, second
## 
## 
## Lade nötiges Paket: grid
## 
## ========================================
## ComplexHeatmap version 2.22.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
## 
## If you use it in published research, please cite either one:
## - Gu, Z. Complex Heatmap Visualization. iMeta 2022.
## - Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
##     genomic data. Bioinformatics 2016.
## 
## 
## The new InteractiveComplexHeatmap package can directly export static 
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
## 
## This message can be suppressed by:
##   suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
## 
## 
## ========================================
## circlize version 0.4.16
## CRAN page: https://cran.r-project.org/package=circlize
## Github page: https://github.com/jokergoo/circlize
## Documentation: https://jokergoo.github.io/circlize_book/book/
## 
## If you use it in published research, please cite:
## Gu, Z. circlize implements and enhances circular visualization
##   in R. Bioinformatics 2014.
## 
## This message can be suppressed by:
##   suppressPackageStartupMessages(library(circlize))
## ========================================
## 
## 
## 
## Attache Paket: 'plotly'
## 
## 
## Das folgende Objekt ist maskiert 'package:ComplexHeatmap':
## 
##     add_heatmap
## 
## 
## Das folgende Objekt ist maskiert 'package:AnnotationDbi':
## 
##     select
## 
## 
## Das folgende Objekt ist maskiert 'package:ggplot2':
## 
##     last_plot
## 
## 
## Das folgende Objekt ist maskiert 'package:IRanges':
## 
##     slice
## 
## 
## Das folgende Objekt ist maskiert 'package:S4Vectors':
## 
##     rename
## 
## 
## Das folgende Objekt ist maskiert 'package:stats':
## 
##     filter
## 
## 
## Das folgende Objekt ist maskiert 'package:graphics':
## 
##     layout
dir.create("results/DESeq", showWarnings = F, recursive = T)

knitr::opts_chunk$set(fig.width = 10, dpi = 300, results = "hold", fig.show = "hold")
## Counts plot

countsToPlot <- c("CD274", "BRAF", "BCL2", "VEGFC", "PDGFC")

## Heatmaps

hm_scale_by_row <- TRUE
hm_max_rows <- 50

hmap_colors <- colorRamp2(c(-1.5, 0, 1.5), c("darkblue", "white", "darkred"))
annotation_colors <- c("nHEM" = "#F8766D",
                       "MCM1G" =  "#7CAE00",
                       "MCMDLN" = "#00BFC4",
                       "MCMDLN_3D" = "#C77CFF")

## Volcano plot

vp_max_labels <- 0
vp_lfc_limit <- NA

Import data

# Read metadata
samples <- read.delim("data/metadata.txt") %>% as_tibble()

# Read featurecounts output
countdata <- read_counts(samples$Files, "data/featurecounts/", sample_IDs = samples$Sample_ID)

# Read GTF file
genedata <- read_genes("data/genes.gtf", countdata$Geneid)
## Automatically detected number of rows to skip:  0 
## List of features in column 9:
## -----------------------------
## gene_biotype 
## gene_id 
## gene_name 
## gene_source 
## transcript_id 
## transcript_name 
## transcript_source 
## tss_id

Prepare data for analysis

countdata.filtered <- filter_counts(countdata)

# Move gene info to separate dataframe and remove from count data
genedata.filtered <- genedata %>% 
  filter(ENSEMBL %in% rownames(countdata.filtered))

samples %<>%
  mutate(Type = str_remove(Sample_ID, "(_[23]D)?_\\d$"),
         Class = factor(Class, c("nHEM", "MCM1G", "MCMDLN", "MCMDLN_3D")))

summary(samples)
##   Sample_ID            Files             Replicate  Annotation       
##  Length:12          Length:12          Min.   :1   Length:12         
##  Class :character   Class :character   1st Qu.:1   Class :character  
##  Mode  :character   Mode  :character   Median :2   Mode  :character  
##                                        Mean   :2                     
##                                        3rd Qu.:3                     
##                                        Max.   :3                     
##      Type            Dimension               Class  
##  Length:12          Length:12          nHEM     :3  
##  Class :character   Class :character   MCM1G    :3  
##  Mode  :character   Mode  :character   MCMDLN   :3  
##                                        MCMDLN_3D:3  
##                                                     
## 

Run DESeq analysis

# Create DESeq dataset
dds <- DESeqDataSetFromMatrix(countData = countdata.filtered, colData = samples, design = ~Class)

# Run DESeq analysis
dds = DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
normalized_counts = counts(dds, normalized = T)

Overview and QC plots

# Size Factors for Normalization
sizeFactors(dds)

# Dispersion plot for QC
plotDispEsts(dds)

# transform data with regularized log for visualization
vsd <- vst(dds, blind = T)

# Principal component analysis
plot_pca(vsd)

# Correlation heatmap
samplecorr <- cor(assay(vsd))
rownames(samplecorr) <- samples$Sample_ID

hmap_annotation <- samples %>%
  dplyr::select(c(Class))

Heatmap(samplecorr, c("darkblue", "khaki1", "darkred"), name = "Correlation", 
        show_column_names = FALSE, row_dend_reorder = FALSE, column_dend_reorder = FALSE,
        top_annotation = HeatmapAnnotation(Class = hmap_annotation$Class,
                                           col = list(Class = annotation_colors),
                                           show_annotation_name = FALSE))
##      nHEM_1      nHEM_2      nHEM_3     MCM1G_1     MCM1G_2     MCM1G_3 
##   0.4848247   0.6124169   0.7109039   1.4445276   1.3167976   1.1891653 
##    MCMDLN_1    MCMDLN_2    MCMDLN_3 MCMDLN_3D_1 MCMDLN_3D_2 MCMDLN_3D_3 
##   0.9962995   1.0757996   1.2974828   1.1285081   1.0287994   1.5023460

Extract results from DESeq object

res.1GvsHEM = results(dds, contrast = c("Class", "MCM1G", "nHEM"), alpha = 0.05)
res.DLNvsHEM = results(dds, contrast = c("Class", "MCMDLN", "nHEM"), alpha = 0.05)
res.DLNvs1G = results(dds, contrast = c("Class", "MCMDLN", "MCM1G"), alpha = 0.05)
res.3Dvs2D = results(dds, contrast = c("Class", "MCMDLN_3D", "MCMDLN"), alpha = 0.05)


plotMA(res.1GvsHEM, ylim=c(-2, 2), main = "MCM1G vs nHEM")
plotMA(res.DLNvsHEM, ylim=c(-2, 2), main = "MCMDLN vs nHEM")
plotMA(res.DLNvs1G, ylim=c(-2, 2), main = "MCMDLN vs MCM1G")
plotMA(res.3Dvs2D, ylim=c(-2, 2), main = "MCMDLN (3D) vs MCMDLN (2D)")


print("MCM1G vs nHEM")
summary(res.1GvsHEM)

print("MCMDLN vs nHEM")
summary(res.DLNvsHEM)

print("MCMDLN vs MCM1G")
summary(res.DLNvs1G)

print("MCMDLN (3D) vs MCMDLN (2D)")
summary(res.3Dvs2D)
## [1] "MCM1G vs nHEM"
## 
## out of 46399 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 9536, 21%
## LFC < 0 (down)     : 11374, 25%
## outliers [1]       : 13, 0.028%
## low counts [2]     : 9896, 21%
## (mean count < 1)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
## 
## [1] "MCMDLN vs nHEM"
## 
## out of 46399 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 10379, 22%
## LFC < 0 (down)     : 11900, 26%
## outliers [1]       : 13, 0.028%
## low counts [2]     : 9896, 21%
## (mean count < 1)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
## 
## [1] "MCMDLN vs MCM1G"
## 
## out of 46399 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 8332, 18%
## LFC < 0 (down)     : 7563, 16%
## outliers [1]       : 13, 0.028%
## low counts [2]     : 10795, 23%
## (mean count < 1)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
## 
## [1] "MCMDLN (3D) vs MCMDLN (2D)"
## 
## out of 46399 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 5219, 11%
## LFC < 0 (down)     : 5392, 12%
## outliers [1]       : 13, 0.028%
## low counts [2]     : 15293, 33%
## (mean count < 3)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

Perform shrinkage

res.1GvsHEM <- lfcShrink(dds = dds, res = res.1GvsHEM, type = "ashr")
## using 'ashr' for LFC shrinkage. If used in published research, please cite:
##     Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2.
##     https://doi.org/10.1093/biostatistics/kxw041
res.DLNvsHEM <- lfcShrink(dds = dds, res = res.DLNvsHEM, type = "ashr")
## using 'ashr' for LFC shrinkage. If used in published research, please cite:
##     Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2.
##     https://doi.org/10.1093/biostatistics/kxw041
res.DLNvs1G <- lfcShrink(dds = dds, res = res.DLNvs1G, type = "ashr")
## using 'ashr' for LFC shrinkage. If used in published research, please cite:
##     Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2.
##     https://doi.org/10.1093/biostatistics/kxw041
res.3Dvs2D <- lfcShrink(dds = dds, res = res.3Dvs2D, type = "ashr")
## using 'ashr' for LFC shrinkage. If used in published research, please cite:
##     Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2.
##     https://doi.org/10.1093/biostatistics/kxw041
plotMA(res.1GvsHEM, ylim=c(-2, 2), main = "MCM1G vs nHEM")
plotMA(res.DLNvsHEM, ylim=c(-2, 2), main = "MCMDLN vs nHEM")
plotMA(res.DLNvs1G, ylim=c(-2, 2), main = "MCMDLN vs MCM1G")
plotMA(res.3Dvs2D, ylim=c(-2, 2), main = "MCMDLN (3D) vs MCMDLN (2D)")

Results

results_all <- list(MCM1G_vs_nHEM = res.1GvsHEM, 
                    MCMDLN_vs_nHEM = res.DLNvsHEM, 
                    MCMDLN_vs_MCM1G = res.DLNvs1G, 
                    MCMDLN_3D_vs_2D = res.3Dvs2D)

results_all %>%
  lapply(function(x) {
    x %>% 
      as.data.frame() %>%
      filter(padj < 0.05 & abs(log2FoldChange) >= 1) %>%
      dplyr::select(log2FoldChange) %>%
      sign() %>%
      mutate(Expression = ifelse(log2FoldChange == 1, "Overexpressed", "Underexpressed")) %>%
      dplyr::select(!log2FoldChange) %>%
      group_by(Expression) %>%
      count()
  }) %>%
  bind_rows(.id = "Comparison") %>%
  pivot_wider(names_from = Expression, values_from = n) %>%
  mutate(Total = Overexpressed + Underexpressed)

Expression levels (all samples)

counts.plot <- lapply(countsToPlot, function(x) {
  plotCounts(dds, genedata.filtered[genedata.filtered$Name == x,]$ENSEMBL,
              intgroup = c("Class"), returnData = TRUE)
  }) %>%
  set_names(countsToPlot) %>%
  lapply(rownames_to_column, var = "Sample_ID") %>%
  bind_rows(.id = "Gene") %>%
  group_by(Gene, Class) %>%
  summarise(mean_count = mean(count), sd_count = sd(count), .groups = "drop") %>%
  mutate(min_error = ifelse(mean_count - sd_count >= 0, mean_count - sd_count, 0), 
           max_error = mean_count + sd_count)

ggplot(counts.plot, aes(x = Gene, y = mean_count, fill = Class)) +
  geom_col(position = "dodge") +
  geom_errorbar(aes(ymin = min_error, ymax = max_error), 
                position = position_dodge(0.9), width = 0.2) +
  labs(x = "", y = "Mean Count", fill = "") +
  theme_pubr() +
  theme(legend.position = "right", legend.text = element_text(size = 12), panel.grid.major.y = element_line(linewidth = 0.5))

MCM1G vs nHEM

res.1GvsHEM.tb <- res.1GvsHEM %>%
  data.frame() %>%
  rownames_to_column(var = "ENSEMBL") %>%
  as_tibble() %>%
  left_join(genedata.filtered, ., by = "ENSEMBL") %>%
  mutate(threshold = padj < 0.05 & abs(log2FoldChange) >= 1)
  
res.1GvsHEM.sig <- res.1GvsHEM.tb %>%
  filter(threshold == TRUE)

pvalues <- res.1GvsHEM.sig$padj
names(pvalues) <- rownames(res.1GvsHEM.sig)

norm.1GvsHEM.sig <- normalized_counts %>%
  data.frame() %>%
  dplyr::select(samples$Sample_ID[samples$Class %in% c("MCM1G", "nHEM")]) %>%
  filter(rownames(.) %in% res.1GvsHEM.sig$ENSEMBL) 

write.csv(res.1GvsHEM.tb %>% dplyr::select(-threshold), 
          file = "results/DESeq/MCM1G_vs_nHEM_unfiltered.csv",
          row.names = FALSE)
write.csv(res.1GvsHEM.sig %>% dplyr::select(-threshold), 
          file = "results/DESeq/MCM1G_vs_nHEM_filtered.csv",
          row.names = FALSE)
plot_heatmap(norm.1GvsHEM.sig, hmap_colors, annotation_colors, 
             genedata.filtered, samples, pvalues, max_rows = hm_max_rows)

volcanoplot(res.1GvsHEM.tb, title = "MCM1G vs. nHEM", max_labels = vp_max_labels, 
            lfc_limit = vp_lfc_limit)

volcanoplot_interactive(res.1GvsHEM.tb, title = "MCM1G vs. nHEM")

MCMDLN vs nHEM

res.DLNvsHEM.tb <- res.DLNvsHEM %>%
  data.frame() %>%
  rownames_to_column(var = "ENSEMBL") %>%
  as_tibble() %>%
  left_join(genedata.filtered, ., by = "ENSEMBL") %>%
  mutate(threshold = padj < 0.05 & abs(log2FoldChange) >= 1)
  
res.DLNvsHEM.sig <- res.DLNvsHEM.tb %>%
  filter(threshold == TRUE)

pvalues <- res.DLNvsHEM.sig$padj
names(pvalues) <- rownames(res.DLNvsHEM.sig)

norm.DLNvsHEM.sig <- normalized_counts %>%
  data.frame() %>%
  dplyr::select(samples$Sample_ID[samples$Class %in% c("MCMDLN", "nHEM")]) %>%
  filter(rownames(.) %in% res.DLNvsHEM.sig$ENSEMBL) 

write.csv(res.DLNvsHEM.tb %>% dplyr::select(-threshold), 
          file = "results/DESeq/MCMDLN_vs_nHEM_unfiltered.csv",
          row.names = FALSE)
write.csv(res.DLNvsHEM.sig %>% dplyr::select(-threshold), 
          file = "results/DESeq/MCMDLN_vs_nHEM_filtered.csv",
          row.names = FALSE)
plot_heatmap(norm.DLNvsHEM.sig, hmap_colors, annotation_colors, 
             genedata.filtered, samples, pvalues, max_rows = hm_max_rows)

volcanoplot(res.DLNvsHEM.tb, title = "MCMDLN vs. nHEM", max_labels = vp_max_labels, 
            lfc_limit = vp_lfc_limit)

volcanoplot_interactive(res.DLNvsHEM.tb, title = "MCMDLN vs. nHEM")

MCMDLN vs MCM1G

res.DLNvs1G.tb <- res.DLNvs1G %>%
  data.frame() %>%
  rownames_to_column(var = "ENSEMBL") %>%
  as_tibble() %>%
  left_join(genedata.filtered, ., by = "ENSEMBL") %>%
  mutate(threshold = padj < 0.05 & abs(log2FoldChange) >= 1)
  
res.DLNvs1G.sig <- res.DLNvs1G.tb %>%
  filter(threshold == TRUE)

pvalues <- res.DLNvs1G.sig$padj
names(pvalues) <- rownames(res.DLNvs1G.sig)

norm.DLNvs1G.sig <- normalized_counts %>%
  data.frame() %>%
  dplyr::select(samples$Sample_ID[samples$Class %in% c("MCMDLN", "MCM1G")]) %>%
  filter(rownames(.) %in% res.DLNvs1G.sig$ENSEMBL) 

write.csv(res.DLNvs1G.tb %>% dplyr::select(-threshold), 
          file = "results/DESeq/MCMDLN_vs_MCM1G_unfiltered.csv",
          row.names = FALSE)
write.csv(res.DLNvs1G.sig %>% dplyr::select(-threshold), 
          file = "results/DESeq/MCMDLN_vs_MCM1G_filtered.csv",
          row.names = FALSE)
plot_heatmap(norm.DLNvs1G.sig, hmap_colors, annotation_colors, 
             genedata.filtered, samples, pvalues, max_rows = hm_max_rows)

volcanoplot(res.DLNvs1G.tb, title = "MCMDLN vs. MCM1G", max_labels = vp_max_labels, 
            lfc_limit = vp_lfc_limit)

volcanoplot_interactive(res.DLNvs1G.tb, title = "MCMDLN vs. MCM1G")

MCMDLN (3D) vs MCMDLN (2D)

res.3Dvs2D.tb <- res.3Dvs2D %>%
  data.frame() %>%
  rownames_to_column(var = "ENSEMBL") %>%
  as_tibble() %>%
  left_join(genedata.filtered, ., by = "ENSEMBL") %>%
  mutate(threshold = padj < 0.05 & abs(log2FoldChange) >= 1)
  
res.3Dvs2D.sig <- res.3Dvs2D.tb %>%
  filter(threshold == TRUE)

pvalues <- res.3Dvs2D.sig$padj
names(pvalues) <- rownames(res.3Dvs2D.sig)

norm.3Dvs2D.sig <- normalized_counts %>%
  data.frame() %>%
  dplyr::select(samples$Sample_ID[samples$Class %in% c("MCMDLN", "MCMDLN_3D")]) %>%
  filter(rownames(.) %in% res.3Dvs2D.sig$ENSEMBL) 

write.csv(res.3Dvs2D.tb %>% dplyr::select(-threshold), 
          file = "results/DESeq/MCMDLN3D_vs_MCMDLN_unfiltered.csv",
          row.names = FALSE)
write.csv(res.3Dvs2D.sig %>% dplyr::select(-threshold), 
          file = "results/DESeq/MCMDLN3D_vs_MCMDLN_filtered.csv",
          row.names = FALSE)
plot_heatmap(norm.3Dvs2D.sig, hmap_colors, annotation_colors, 
             genedata.filtered, samples, pvalues, max_rows = hm_max_rows)

volcanoplot(res.3Dvs2D.tb, title = "MCMDLN 2D vs. 3D", max_labels = vp_max_labels, 
            lfc_limit = vp_lfc_limit)

volcanoplot_interactive(res.3Dvs2D.tb, title = "MCMDLN 2D vs. 3D")